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FOREWORD 


This  report  presents  a  summary  of  work  performed  for  the  Air  Force  Flight 
Dynamics  Laboratory,  Air  Force  Systems  Command,  funded  under  Contract 
F33615-72-C-1429.  The  work  was  performed  by  the  Aerodynamics  Research 
Department  of  the  Northrop  Corporation,  Aircraft  Division,  Hawthorne,  California. 

This  report  is  divided  into  two  volumes.  Volume  I  discusses  the  theory  and 
the  application,  and  presents  comparisons  of  the  numerical  results  with  experi¬ 
mental  data.  Volume  TT  discusses  the  details  of  the  computer  programs  and  how 
to  operate  them.  For  internal  control  purposes,  this  report  has  been  assigned 
the  Northrop  report  number,  NOR  72-87. 

The  work  reported  herein  was  begun  on  December  1,  1971,  and  the  final,  opera 
tional  computer  programs  were  delivered  and  demonstrated  in  the  period  of  -Tune 
14-16,  1972.  In  addition  to  the  two  authors,  Mr.  Joe  Der,  Jr. ,  also  of  the  Aero¬ 
dynamics  Research  Department,  made  significant  contributions  to  this  work 
through  his  knowledge  of  computing  and  computing  machinery.  Mr.  Henry  Ziegler 
did  a  fine  job  in  the  preparation  of  the  final  report,  and  Miss  Alberta  Hansen  made 
all  the  line  drawings. 

Mr.  A.  B.  Lewis  was  the  Air  Force  Project  Engineer.  He  also  made  very 
important  contributions  to  the  program,  and  much  of  the  credit  for  the  back-to  - 
back  operation  of  the  two  programs  must  go  to  him  for  his  suggestions  and  gui¬ 
dance. 


This  technical  report  was  submitted  by  the  authors  in  July  1972,  and  has  been 
reviewed  and  approved. 
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SECTION  I 


INTRODUCTION 

Tbe  inviscid  flow  field  programs  contained  within  the  combined  viscid-inviscid 
flow  field  program  described  in  Reference  1  have  been  separated  from  the  overall 
program,  revised  extensively  and  improved  in  many  respects.  The  logic  of  these 
programs  has  been  greatly  simplified;  unused  variables  and  unnecessary  blocks 
of  logic  have  been  removed,  and  many  new  developments  added.  As  a  result,  two 
computer  programs  are  now  available  for  calculating  the  inviscid  supersonic/hyper¬ 
sonic  flow  over  three-dimensional  bodies  at  angles  of  attack  with  a  high  degree  of 
accuracy  and  for  a  modest  expenditure  of  computer  time. 

The  resulting  new  programs  are  the  Initial  Value  Surface  Program  and  the 
Three-Dimensional  Method  of  Characteristics  Program.  Both  are  based  on  the 
use  of  an  ideal  gas  with  a  constant  ratio  of  specific  heats.  Both  programs  have 
been  designed  to  require  a  minimum  of  input  and  to  be  extremely  flexible  to  use. 

All  development  decisions  were  based  on  the  use  of  both  programs  by  the  Air  Force 
Flight  Dynamics  Laboratory  on  the  Wright-Patterson  Air  Force  Base  CDC  6600 
computer;  use  on  other  machines  may  result  in  some  compromises. 

The  bodies  which  may  be  .reated  must  have  spherical  noses  and  no  slope  dis¬ 
continuities,  but  may  have  upper  and  lower  flat  sections  (for  example,  a  slab  delta 
wing).  Any  angle  of  attack  for  which  the  subsonic  region  does  not  extend  beyond 
the  spherical  nose  can  be  handled. 

Volume  I  discusses  the  theoretical  founds  and  the  numerical  methods 
while  Volume  II  describes  the  computer  programs  and  how  to  use  them.  The 
original  methods  were  given  in  Reference  1;  a  complete,  up-to-date  discussion  of 
the  basic  formulation  and  the  numerical  methods  is  presented  here. 
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SECTION  II 


BODY  DESCRIPTION 

Figure  1  shows  a  complex  body  shape  having  a  spherical  nose  and  upper  and 
lower  flats.  In  a  right-handed  coordinate  system,  the  Y-axis  is  aligned  with  the 
body  axis.  The  Z-axis  is  up  and  the  X-axis  spanwise.  The  nose  of  the  body  need 
not  coincide  with  the  origin  of  the  coordinate  system. 

Three  lines  serve  to  define  this  body  shape: 

(1)  The  outer  limit  of  the  upper  flat 

(2)  The  maximum  width  line 

(3)  The  outer  limit  of  the  lower  flat. 

In  cross  section,  an  ellipse  is  fitted  between  line  (1)  and  (2),  and  another 
ellipse  between  (2)  and  (3) .  For  the  case  of  a  blunted  cone,  (or  ellipse)  the  outer 
limits  of  the  upper  and  lower  flats  lie  in  the  plane  of  symmetry  X  *  0. 

Each  line  is  described  by  two  functions  of  Y. 

H,  -  T?Y  +  Q,  +  Sfi,yR,Y*+  S,Y  +  T, 

X,  -  P«Y  +  04*  S,Y+T4 

za-  Pj  y  +  q2  +  ssay/?aYl+  s,Y+r, 

Xj-  P,  Y  +  Qs  SS3y  R,Y*+  S5V+  r/ 
z5-  P,r  +Qj  +  SSjV  P*Y*+  S,Y+T, 

X,=  P,Y  +  a*  +  S6*y  ff,Y*  +  S»Y  ♦  T* 

Here  P,  Q,  R,  S  and  T  are  coefficients  of  the  general  conic  sections  and  SG  is  the 
sign  of  the  square  root.  Methods  for  evaluation  of  these  coefficients  are  discussed 
in  Appendix  I. 

No  slope  discontinuities  are  permitted  in  the  present  program.  Since  each 
line  may  be  described  by  several  segments,  this  slope  smoothness  requirement 
translates  into: 


mi 

everywhere  continuous 

m 

everywhere  continuous 

m 

everywhere  continuous 

However, 


dX, 

dY 


dlt 

dY 


and 


are  not  required  to  be  continuous. 


Two  examples  of  this  body  description  method  are  given  in  Appendix  II. 


SECTION  HI 

THE  INITIAL  VALUE  SURFACE 

The  function  of  the  Initial  Value  Surface  (IVS)  Program  is  to  provide  sufficient 
data  to  start  the  Three-Dimensional  Method  of  Characteristics  Program  calculation 
of  the  supersonic  flow  field.  For  reasons  of  economy,  this  IVS  should  be  as  far 
downstream  as  possible  in  order  to  compute  the  largest  possible  portion  of  the 
flow  field  by  means  of  a  rotationally  symmetric  method.  This,  of  course,  implies 
a  minimization  of  the  total  computing  time  required  for  the  total  solution. 

The  present  IVS  program  is  a  greatly  modified  version  of  the  old  IVS  program 
from  Reference  1.  Figure  2  compares  the  old  and  new  approaches. 

Since  a  sphere  has  no  preferential  direction,  the  solution  for  a  sphere  at  angle 
of  attack  can  be  obtained  by  rotating  the  solution  for  a  sphere  at  zero  angle  of 
attack.  The  old  IVS  was  obtained  by  calculating  a  zero  angle  of  attack  blunt  body 
solution,  forming  an  initial  value  line  (IVL),  calculating  the  supersonic  flow  field 
up  to  the  outgoing  characteristic  from  the  foot  of  the  IVL,  rotating  this  limiting 
characteristic  around  the  free  stream  velocity  vector,  and  pitching  over  to  the 
correct  angle  of  attack. 

The  new  IVS  program  works  in  a  similar  manner  up  through  the  generation  of 
the  IVL.  After  the  IVL  is  generated,  the  Rotationally  Symmetric  Method  of  Char¬ 
acteristics  solution  is  carried  sufficiently  far  around  the  sphere  so  that  data  along 
a  series  of  lines  (one  line  for  each  meridian  plane)  can  be  determined  by  interpo¬ 
lation;  these  lines  are  then  rotated  around  the  freestream  vector  and  the  resulting 
plane  pitched  over  by  the  angle  of  attack.  The  resulting  plane  is  now  perpendicular 
to  the  body  axis  and  situated  at  the  sphere-body  juncture  (see  Figure  2).  In  theory, 
the  shaded  region  shown  in  Figure  2  could  also  be  calculated  by  rotationally  symmet¬ 
ric  means  but  the  resulting  surface  geometry  is  too  complex  to  permit  general  use. 
Therefore,  the  constant  body  station  IVS  represents  a  good  compromise,  and  is  a 
great  improvement  over  the  previous  program. 

The  IVS  program  consists  of  four  separate  programs,  all  connected  together 
through  an  OVERLAY  sequence.  These  programs  are: 
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FIGURE  2.  COMPARISON  BETWEEN  OLD  AND  NEW  INITIAL  VALUE  SURFACES 


(1)  Blunt  Body  Program 

(2)  Initial  Value  Line  Interpolation  Program 

(3)  Rotationally  Symmetric  Method  of  Characteristics  Program 

(4)  Initial  Value  Surface  Interpolation  Program 

The  Blunt  Body  Program  determines  the  subsonic  through  slightly  supersonic 
flow  over  the  blunted  nose.  The  Initial  Value  Line  Interpolation  Program  inter¬ 
polates  in  the  Blunt  Body  Program  data  to  develop  a  start  line  for  the  Rotationally 
Symmetric  Method  of  Characteristics  Program.  This  latter  program  calculates 
the  required  remaining  supersonic  flow  field  over  the  nose.  The  Initial  Value 
Surface  Interpolation  Program  interpolates  in  the  RSMoC-generated  data  to  deter¬ 
mine  properties  on  a  plane  normal  to  the  body  axis,  given  the  number  of  data  rings 
and  the  number  of  meridian  planes.  Each  of  these  programs  is  discussed  below. 

1.  BLUNT  BODY  PROGRAM 

The  theoretical  foundations  of  this  program  were  discussed  in  Reference  1.  No 
significant  logic  changes  were  made  to  the  program,  and  so  the  following  discussion 
is  essentially  a  summary  of  the  approach. 

The  flow  field  behind  the  detached  bow  shock  wave  of  the  spherical  nose  is  cal¬ 
culated  by  an  extension  of  Van  Dyke's  method  (Reference  2),  and  is  an  indirect  one 
(i.  e. ,  the  calculation  proceeds  from  a  specified  shock  to  an  unknown  body) .  The 
shape  of  this  unknown  body  is  compared  with  that  of  the  desired  body,  and  appropri¬ 
ate  changes  are  made  in  the  shock  wave  shape  until  the  difference  between  calcula¬ 
ted  and  desired  bodies  is  acceptable. 

The  calculation  is  carried  out  in  a  curvilinear  orthogonal  coordinate  system 
based  on  the  assumed  shock  wave.  The  shock  wave  is  taken  as  a  general  conic 

Ya»  2RX-BX*  (l) 

where  X  and  Y  are  the  usual  Cartesian  coordinates,  R  is  the  radius  of  curvature  of 
the  shock  wave  at  the  origin,  and  B  is  the  "bluntness."  The  shock  is  hyperbolic  for 
B<0,  parabolic  for  B=0,  and  elliptic  for  B>0.  By  comparing  Equation  (1)  with  the 
standard  form  of  the  conic  equation,  it  can  be  seen  that  for  B>0,  B  =b2/a2  ;  where 
b  is  the  semi-axis  in  the  Y  direction,  and  a  the  semi-axis  in  the  X  direction.  There¬ 
fore,  B  =  1  corresponds  to  a  circle  or  a  sphere. 
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a 


A  curvilinear  orthogonal  coordinate  system,  §  ,  rj ,  based  on  the  shock  wave 
shape  of  Equation  (1)  then  can  be  erected.  The  relationship  between  the  Cartesian 
and  curvilinear  coordinates  can  be  shown  to  be 


«*  ~  (2) 

and 

•£-  Iflt?  <3> 

where  the  shock  wave  corresponds  to  =  1.  Figure  3  illustrates  this  coordinate 
system. 


Y 


x 

FIGURE  3.  CURVILINEAR  COORDINATE  SYSTEM 

To  circumvent  certain  singularities  in  the  governing  gas  system  of  dynamics 
equations,  the  first  §  used  by  the  finite-difference  solution  is  taken  to  be  Af/2; 
for  succeeding  values,  §n  +  A  f  .  Note  that  the  body  will  not  corres¬ 

pond  to  a  constant  surface. 

By  introducing  a  Stokes*  stream  function,  V , 


A 

3 


i 

3 

■> 


P/p *  =  f(Y) 


(4) 


the  governing  gas  dynamics  equations  can  be  shown  to  be 


(6) 


Note  that,  while  P,  is  a  function  of  Vf  ,  ,  fa,  ,  and  Pj  ,  %n  ls  a 

function  of  ,  Vfy ,  P|  ,  and  P^ .  Therefore  Equation  (5)  must  be 

evaluated  to  provide  P^  for  Equation  (6). 
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In  Equations  (5)  and  (6),  S'/R  denotes  jjy  (since  S  is  only  a  function  of  ¥ 
in  the  fluid  mechanics  sense).  For  an  ideal  gas,  the  entropy  derivatives  reduce  to: 


(Z(6/R)\  _  _v_  1 

V  ap  Jp  v-i  p 

(7) 

/a(s/R)\  ^  i  i 
\  qp  /p  v-i  p 

(8) 

and 

d(S/ R)  i  f 

dv  v-i  T 

0) 

where 

P/p'  -  f. 

The  computation  scheme  is  as  follows: 

1. 

On  the  shock,  determine  all  properties  at  a  series  of  evenly  spaced  points 

2. 

Determine  the  §  derivatives. 

3. 

Determine  tb"  Pq  and  values  from  Equations  (5)  and  (6) . 

4. 

Integrate  into  the  field  one  Arf  step  by  using  the  following: 

P2  -  P/  - 

(10) 

(11) 

and 

(12) 

5.  Knowing  P  and  at  the  new  jp  coordinate,  evaluate  all  the  necessary 
properties. 

6.  By  numerical  means,  determine  the  §  derivatives  of  P,  ¥ ,  V ,  {^  , 
and  ¥§  at  the  new  coordinate  line. 

7.  Calculate  P^  and  at  the  new  coordinate  line,  using  Equations  (5) 

and  (6). 
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8.  Determine  the  appropriate  physical  properties  at  the  new  coordinate  line, 
by  an  "external  iteration"  on  the  pressure  through  the  following  equation: 

Pa  -  Pf  (Pn  +  Pija')  (13) 

Note  that  is  used  only  to  evaluate  the  flow  field  properties  at  the  new  point, 

and  not  to  reevaluate  Pjj  or  V^n  • 

9.  Thus  the  properties  and  derivatives  at  the  new  coordinate  line  are  known, 
and  this  line  then  may  be  used  as  a  base  from  which  another  new  coordi¬ 
nate  line  can  be  calculated. 

10.  This  scheme  is  continued  until  negative  values  of  the  stream  function  are 
obtained.  The  body  location  is  found  by  extrapolating  data  from  the  last 
three  surfaces  having  V>0  to  determine  where  0. 

Using  the  marching  technique  previously  outlined  in  steps  1-10,  the  flow 
field  from  the  shock  to  the  body  can  be  determined.  The  actual  calculation  of  the 
body,  however,  requires  special  techniques.  The  marching  technique  is  used 
until  values  of  V  less  than  zero  appear.  It  has  been  found  that  ,  and 

y\*i  all  change  drastically  inside  the  body,  and  cannot  be  used  to  determine  the 
body  location  by  interpolation.  Therefore,  data  from  the  last  three  =  constant 
lines  are  used  to  extrapolate  to  the  Y-0  point  to  determine  the  body.  For  each 
V  «  constant  line,  7?  body,  P,  Yj  and  also  are  extrapolated  as  a  func- 


a  normal  shock,  all  properties  at  the  body  surface  can  be  determined. 


Figure  4  illustrates  a  complete  solution  for  M®  =  8. 8  in  air.  Note  that  each 
step  requires  the  dropping  of  six  end  points  from  the  analysis.  This  is  a 
result  of  using  the  finite-difference  differentiation  technique.  To  determine  Pj  , 
three  points  on  either  side  must  be  used.  Hence,  the  last  P|  that  can  be  deter¬ 
mined  is  the  fourth  from  the  end  of  the  If  =  constant  line.  After  the  first  forward 
integration  has  occurred,  P|  on  the  new  line  must  be  determined;  again,  the  last 
three  points  are  required  to  determine  the  value  for  the  fourth  from  the  end.  Thus 
a  total  of  six  points  must  be  dropped  for  every  step  taken  into  the  flow  field. 

The  determination  of  the  body  shape  is  based  on  the  assumption  that  a  conic 
shock  will  yield  a  conic  body.  Hence  the  body  equation  can  be  written  as: 
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FIGURE  4.  COMPLETE  FLOW  FIELD  FOR  A  SPHERE 


(14) 


Yg  =  2/?g(Ag-4) 

where  R„  =  radius  of  curvature  of  the  body  at  the  origin,  referenced  to  R 
of  the  shock  wave 

Bg  =  body  bluntness 

A  =  shock  stand-off  distance,  referenced  to  R  of  the  shock  wave 

To  determine  the  best  values  for  Rg,  Bg,  and  A  for  the  N  body  points  deter¬ 
mined  by  the  program,  the  following  set  of  simultaneous  equations  (which  represent 
the  least-square  solution  to  the  body  description  equation  with  an  off-set  axis)  must 
be  solved  for  <X.,  (&,  and  V: 

on  +  jSExi  +  *£Xi4  *  EVi*  as) 

otExi  +  $  E  x  *  +  *EXis  -  EXvVi2  (i6) 

oc£Xi*  ♦  (&£xi3  v£Xi4  -  EXiaY\*  a?) 

where  E  denotes  summation  from  i  =  1  to  N,  and  X.,  Y.  are  the  Cartesian  coordi¬ 
nates  of  the  body  points. 

The  body  parameters  then  are  determined  from  the  following  functions  of  <X, 

/3  ,  and  V: 

<x  --M2R a  +BaA)  (18) 

*»  2(R g  +  BqA)  (19) 

y  -  J&B  (20) 

The  summations  indicated  in  Equations  (15),  (16),  and  (17)  include  only  subsonic 
body  points.  It  has  been  found  that  the  supersonic  body  points  tend  to  deviate  from 
a  conic  section  rapidly,  and  that  a  more  consistent  body  description  can  be  obtained 
by  excluding  them.  Similarly,  it  has  been  found  that  if  the  body  point  data  are  ob¬ 
tained  from  one  or  more  AtJ  steps,  a  poor  body  description  results.  This  diffi¬ 
culty  can  be  eliminated  by  providing  a  means  for  encompassing  all  the  body  points 
with  one  AtJ  step.  Thfs  is  not  to  say  that  by  arranging  for  a  one- ATJ -step  body  the 
accuracy  is  increased.  What  is  obtained  is  a  smoother,  more  uniform  variation 
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of  Bg  with  B,  which  permits  the  use  of  an  automated  itt, ration. 

Let  the  step  from  the  shock  to  the  body  stagnation  point  be  ATf  stag’ 
step  from  the  body  stagnation  point  to  the  body  sci.ic  point  be  Al^g,  and  the 
ratio  of  AtJ's  for  successive  steps  be  G  (a  constant) .  , Elementary  series  summa¬ 
tion  analysis  then  shows  chat  if  n  is  the  number  of  steps  to  the  body1  the  following 
equation  can  be  used  to  determine  G: 


Knowing  and  Aljg  and  specifying  n,  Equation  (21)  may  be  solved  by 

Newton's  method  for  G.  The  initial  value  of  is  determined  by: 

A7?o  -  ^  (22, 

In  practice,  At}gtag  and  48JB  are  unknown  initially.  A  first  problem  is  cal¬ 
culated  with  G  =  1 . 0  and  =  AT?  =  0.02  to  determine  where  the  body  lies.  Using 
data  from  this  run,  and  extrapolating  for  the  body  sonic  point  if  necessary,  a  first 
approximation  for  A1?stag  and  At?B  can  be  obtained.  With  each  successive  cal¬ 
culation,  new  values  of  the  increments  are  calculated  so  as  to  provide  the  best 
possible  fit  to  the  body. 

Certain  safety  factors  are  built  into  the  computer  program  to  ensure  that  all 
body  points  lie  in  the  same  ATI  step.  The  AtJ  step  to  the  stagnation  point  is 
reduced  by  two  percent,  and  the  distance  to  the  sonic  point  increased  by  30  percent. 
Hence: 


Atfshag  ( 1.0  -  Tfsfag  )  0.93 


(23) 


and 

=  “  AVsiaq  ”  *?somc  ^  i  i0  (24) 


These  safety  factors  are  not  optimum;  they  merely  work. 


’Note  that  if  n  is  the  number  of  steps  from  the  shock  to  the  body,  there  will  be  n+1 
constant  lines  plus  the  body,  or  n+2  series  of  data. 
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The  overall  iteration  scheme  for  determining  the  shock  wave  corresponding 
to  a  given  body  is  based  on  the  fact  that  Bg  is  a  monotonically  increasing  function 
of  B,  as  shown  in  Figure  5. 


FIGURE  5.  BODY  BLUNTNESS  AS  FUNCTION  OF  SHOCK  BLUNTNESS 

The  iteration  process  takes  the  first  computed  body  bluntness  and  compares 
it  with  the  desired  bluntness.  If  the  calculated  bluntness  is  not  sufficiently  close 
to  the  desired  value  the  shock  bluntness  is  adjusted,  for  the  first  pass  only,  by 
the  following: 

B2  =  B,  (  )  (25) 

For  succeeding  calculations  the  present  and  past  (B,  Bg)  pairs  are  retained 
and  the  new  B  determined  by: 

B*  "  bb1”  Bflj  ^Bb wanked”  BB*)  +  B2  (26) 

The  value  of  Bg,  as  determined  by  a  least  squares  curve  fit,  is  a  weak  function  of 
the  number  of  points  used  in  the  curve  fit.  During  the  iteration  process,  successive 
body  descriptions  differ  in  the  number  of  body  points  obtained  until  convergence  is 
near.  As  a  result,  there  are  a  number  of  essentially  parallel  curves  in  Figure  5, 
one  for  each  different  number  of  body  points.  In  passing  from  one  curve  to  another, 
the  quantity  (B2  -B1)/(Bg  -Bg  ),  the  calculated  local  slope,  may  be  negative. 


In  this  case,  the  iteration  process  governed  by  Equation  (26)  no  longer  provides  a 
rational  answer.  When  the  quantity  AB/ABg  is  negative,  the  calculation  procedure 
then  merely  adds  an  increment  to  the  existing  value  of  B  equivalent  to  that  used  in 
the  last  pass: 

J3$  88  Bj  +  |  Bf  “  B2  |  (27) 

This  permits  the  iteration  procese  to  continue,  and  provides  a  chance  for  the  nat¬ 
ural  convergence  process  to  come  into  play. 

As  noted  above,  a  first  problem  is  run  to  determine  the  location  of  the  body. 
The  B  for  this  run  is  calculated  by  an  approximate  curve  fit  to  the  data  given  in 
Reference  2. 

B  *  0.5  -  0.S59  (Moo-lOy 119  (28) 

After  the  first  problem  is  run  with  a  G  =  1.0,  a  second  is  run  with  a  calculated  G 
but  the  same  B.  Thereafter,  the  iteration  process  discussed  above  takes  over. 

The  blunt  body  computer  program  is  subject  to  certain  limitations  which  are 
described  below: 

1.  The  number  of  shock  points  is  limited  to  100. 

2.  The  number  of  steps  to  the  body  is  limited  to  less  than  10  but  more  than  2. 

3.  A  limit  of  9  overall  calculations  is  imposed  on  each  run,  due  to  instability 
of  the  method.  For  normal  use,  6  is  recommended. 

4.  For  shock  waves  with  B>0,  I2m*x  *  [fo~k>4t]  must  be  less  than 

At  this  point  the  shock  wave  has  degenerated 
into  a  Mach  wave,  and  further  shock  angle  decay  is  forbidden.  If  this 
occurs,  the  program  automatically  reduces  the  value  of  (normally 
equal  to  0.02)  to  accommodate  the  desired  number  of  points.  In  such 
eases  it  is  recommended  that  the  number  of  shock  points  be  reduced  to 
maintain  a|*<0.02. 

5.  In  the  supersonic  region  of  the  flow  field,  the  denominator  of  Equation  (5), 


eventually  goes  down  to  zero.  If  this  occurs  within  the  region  of  calcu¬ 
lation,  the  data  beyond  this  point  are  meaningless  and  are,  therefore, 
dropped.  This  will  cause  the  program  occasionally  to  drop  more  than  1 
6  points  on  a  given  line.  1  , 

i 

2.  INITIAL  VALUE  LINE  INTERPOLATION  PROGRAM  «' 


Upon  completion  of  the  flow  field  evaluation  portion  of  the  Blunt  Body  Program 
all  of  the  flow  field  data  points  are  available  on  TAPE  4  for  use  ,in  determining  the 
IVL.  Using  these  data,  this  program  will  'determine  data  on  a  straight  line  from 
the  body  to  the  shock  such  that  the  characteristic  directions  at  every  point  lie  on 
the  downstream  side  of  the  IVL.  1  , 


One  significant  change  was  worked  into  the  IVL  program.  In  the  earlier  pro¬ 
gram,  a  value  of  1.05  was  specified  internally  for  ,  the  Mach  number  on  the 

body  at  the  foot  of  the  IVL.  This  quantity  is  now  an  input  variable'.  A  value  of 
1.05  is  recommended  for  all  Mach  numbers  greater  than  2. 5,  and  a  value  of  1. 10 
for  Mach  numbers  less  than  2. 5.  The  large  Mbody  at  the  lower  Mach  numbers 
is  due  to  the  poor  definition  of  the  flow  properties  by  the  Blunt  Body  Program  in 
the  very  low  transonic  range.  i 


The  sonic  line  on  a  rotationally  symmetric  blunt  body  takes  one  of  two  basic 
shapes,  depending  on  whether  the  freestream  Mach  number  is  greater  or  less 
than  approximately  5.0.  These  shapes  are  illustrated  in  Figure  6. 


FIGURE  6.  SONIC  LINE  SHAPES 
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The  IVL  program  examines  the  first  supersonic  point  on  each  of  the  "shock - 
like"  coordinates  and  determines  which  of  these  will  minimize  the  slope  (which  is 
negative)  of  a  line  between  each  of  these  points  and  the  center  of  the  body.  The 
program  then  determines  the  approximate  slope  of  the  sonic  line  at  the  body  by 
using  the  location  of  the  body  sonic  point  and  the  sonic  point  in  the  field  nearest  to 
the  body.  If  this  latter  slope  is  more  positive  than  the  previously  discussed  slope, 
a  high  Mach  number  configuration  exists;  if  the  latter  slope  is  the  more  negative 
of  the  two,  a  low  Mach  number  configuration  exists. 

For  the  high  Mach  number  case,  the  line  from  the  center  of  the  body  passing 
through  the  point  on  the  body  where  M  =  is  chosen  as  the  IVL;  the  data 

interpolation  takes  place  along  this  line. 

For  the  low  Mach  number  case,  the  sonic  line  is  highly  curved  and  is  actually 
made  up  of  the  outer  portion  of  the  right-running  limiting  characteristic  and  the 
inner  portion  of  the  left-running  limiting  characteristic.  In  this  situation,  if  the 
high  Mach  number  method  of  determining  the  IVL  were  to  be  used,  the  IVL  would 
pass  through  the  subsonic  region.  Thus,  a  data  line  with  a  different  slope  and  a 
different  intercept  point  on  the  body  is  needed.  This  slope  is  obtained  by  taking 
the  slop'e  of  the  innermost  segment  of  the  sonic  line  (the  calculation  of  which  has 
been  discussed  above)  and  rotating  it  five  degrees  clockwise. 

i 

A  table  of  body  point  coordinates  and  associated  left-running  characteristic 
angles,  <5  +  M  (where  &  is  the  local  flow  angle  and  jjl  the  Mach  angle),  is  prepared 
to  go  from  the  body  point  data  and  an  interpolation  performed  to  determine  the 
location  of  the  body  point  whose  left-running  characteristic  direction  corresponds 
to  the  IVL  slope.  Having  thus  defined  a  point  on  the  IVL  and  its  slope,  the  inter¬ 
polation  for  the  data  points  throughout  the  flow  field  can  be  carried  out.  However, 
the  body  point  determined  by  this  method  is  not  used  as  a  point  on  the  IVL.  Due  to 
the  drastic  increase  in  error  in  calculating  the  local  flow  angles  near  the  body  for 
low  Mach  numbers,  it  has  been  found  that  better  downstream  calculations  are 
obtained  if  the  solution  goes  directly  from  the  first  field  point  to  the  body,  thus 
ignoring  the  IVL  body  point. 

As  a  final  step,  the  freestream  quantities,  the  body  geometric  quantities  (Rg, 
Bg,  4  ),  the  location  and  flow  properties  of  each  point  on  the  IVL,  and  the  mass- 
entropy  table  are  printed  and  stored  on  tape  for  the  Rotationally  Symmetric 
Method  of  Characteristics  Program. 
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The  mass-entropy  table  is  a  listing  of  shock  ordinate  values  and  the  associated 
total  pressure  ratio  (or  downstream  entropy),  starting  at  the  stagnation  line  and 
including  as  the  last  point  the  shock  point  of  the  IVL.  This  table  is  used  with  a 
characteristics  program  to  ensure  an  exact  determination  of  the  local  total  pressure 
(or  entropy)  in  rotational  method  of  characteristics  procedures  (Reference  3) . 

/igure  7  shows  the  body  location  and  s.-e,  and  the  interpolated  sonic  line  and 
initial  value  line  for  the  sphere  solution  shown  earlier  in  Figure  4.  In  the  final 
output,  after  convergence,  all  the  blunt  body  data  are  normalized  by  R^g^;  there¬ 
fore  the  body  radius  is  here  found  to  be  1.0. 

3.  ROTATIONALLY  SYMMETRIC  METHOD  OF  CHARACTERISTICS  PROGRAM 

Using  the  initial  value  line  as  starting  data  the  Rotationally  Symmetric  Method 
of  Characteristics  (RSMoC)  program  calculates  the  supersonic  flow  over  the  sphere- 
cone  body.  These  supersonic  calculations  are  carried  out  along  characteristics 
starting  at  the  shock  wave  and  traveling  in  toward  the  body  as  well  as  downstream. 
The  calculation  continues  until  there  are  at  least  two  data  points  on  each  character¬ 
istic  downstream  of  the  limit  line,  and  one  shock  point  is  downstream  of  the  limit 
line.  The  limit  line  is  a  line  whose  slope  is  tan  (  7T-  Ot )  and  which  passes  a  dis¬ 
tance  of  (1-sin upstream  of  the  center  of  the  sphere  (see  following  section  on 
the  IVS). 

The  present  RSMoC  program  represents  a  considerable  modification  of  the 
previous  RSMoC  program  described  in  Reference  1.  Some  of  the  major  changes 
are: 

(1)  A  mesh  control  scheme  has  been  added  to  permit  the  calculation  of  charac¬ 
teristics  emanating  from  shock  points  and  to  control  the  step  size  achieved 
on  the  body. 

(2)  A  shock  point  routine  using  shock  segments  as  prescribed  by  the  mesh 
control  program. 

(3)  Inclusion  of  the  limit  line  to  terminate  the  construction  of  characteristics. 

(4)  A  body  point  routine. 

The  basic  numerical  methods  are  discussed  below: 


a.  Field  Point 


The  field  point  solution  in  the  rotationally  symmetric  method  of  characteristics 
procedure  determines  the  properties  at  the  intersection  of  two  characteristics  of 
opposite  families,  using  the  mass-entropy  method  (Reference  3) .  A  schematic  of 
two  characteristics  is  shown  in  Figure  8. 


FIGURE  8.  FIELD  POINT  CALCULATION 


Points  A  and  B  in  Figure  8  are  known  points  in  the  flow  field,  and  point  C  is  to  be 
determined.  The  characteristics  connecting  the  "base  points"  (points  A  and  B) 
with  the  unknown  field  point  are  approximated  by  straight  lines  having  a  slope 
determined  by  averaging  the  characteristic  direction  at  the  known  base  point  with 
the  best  available  value  at  the  unknown  field  point.  To  start  the  process,  the 
properties  at  C  are  assumed  to  be  an  average  of  those  at  A  and  B.  Let  point  B  be 
the  point  situated  on  the  right-running  characteristic  which  lies  in  the  b-ji  direc¬ 
tion  (  6  is  the  local  flow  direction  and  jx  is  the  local  Mach  angle) ,  and  let  A  be 
the  point  on  the  left- running  characteristic  which  lies  in  the  b+jx  direction. 

The  characteristic  passing  through  B  can  be  approximated  by 

7^-  -  tan(<Wlac)  =  Tb  (29) 
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where 


and  the  characteristic  passing  through  A  can  be  approximated  by 


Yc-Va 


tan 


(30) 


The  intersection  then  is  defined  by  eliminating  Yc  from  Equations  (29)  and  (30): 

V  _  YA-YB-(XATA-XBTB)  ,31, 

c  re-TA 

Equation  (31),  however,  cannot  provide  accurate  answers  when  either  or  Tg, 
or  both,  become  very  large.  In  these  cases,  Equations  (29)  and  (30)  can  be  re¬ 
written  to  use  the  cotangent  function  instead  of  the  large  values  of  the  tangent,  and 
the  equivalent  of  Equation  (31)  can  be  derived. 

The  compatibility  equation  (Reference  4)  can  be  written  as: 


d  tri?  _  d6  +  am  6  dL 
y  svijulCosjo.  cos  jjl  V 


0 


on  6  ~JX  (32) 


d  In' P  +  d  6  +  5in6  dL 

V  sm/jloosm  co$m  Y 


on  6  ♦  M  (33) 


where  dLj  is  a  distance  along  the  (  6  -M)  characteristic  and  dLg  is  a  distance 
along  the  (  &  +jx)  characteristic.  Equations  (32)  and  (33)  can  then  be  written  in 
finite-difference  form  as: 


and 


y 


6c  “6b 

- —Be"  ~ 

sinjn  cosjul 


+ 


?IC 

sin  6 
—  Z-%2 
cosm 


M  CPc/Pa)  .  &c ~ aLac 

y  swimACcoSmAC  cosjd.^ 


0  (34) 


0  (35) 
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where  ALgC  ■  y  (*B  *c T* 

4L*c  -  V(Xa-X^*v(Ya-Yc)* 

Solving  Equations  (34)  and  (35)  for  bc , 

4  “  i  |^Pa  -tnP8  +  2,  +  _^j 


gnrS  ^Lge  sin  6  <il-Ae 
cosjx*^  '  cosli*?* 


l/[— 4, 

J/  Lsw2ab 

* 


sw2iT 


Similarly,  solving  for  JKn  Pc , 


i«Pc 


Pc  =  [gtniM^MPa  +  sin2jui*ci«pAj  +  dA  -  6»  - 

» 


;j^+  sm2ja4C]J- 


The  actual  computation  sequence  is  as  follows: 

1.  The  properties  at  point  C  are  set  equal  to  averages  of  those  at  points 
A  and  B. 

2.  The  intersection  of  the  approximations  to  the  characteristics  is  obtained 
from  Equation  (31)  or  its  equivalent. 

3.  6c  is  obtained  from  Equation  (36). 

4.  Both  sides  of  Equation  (37)  are  evaluated,  and  the  difference  between 
left  and  right  sides  is  taken  as  a  residual.  If  the  absolute  value  of  this 
residual  is  less  than  a  prescribed  amount,  the  solution  is  acceptable. 
Otherwise  the  iteration  is  continued. 

5.  If  the  solution  is  not  acceptable,  the  new  value  of  P  is  obtained  from 

v 

Equation  (37). 


The  program  retains  corresponding  values  of  Pp  and  its  associated  residual 
for  the  most  recent  positive  and  negative  residual;  a  linear  interpolation  is  car¬ 
ried  out  for  the  value  of  Pp  corresponding  to  a  zero  residual.  This  process  has 
been  found  to  speed  up  the  convergence  process. 


b.  Shock  Point 


The  calculation  of  the  extension  to  the  bow  shock  wave  uses  the  shock  point  and 
the  interpolated  field  point  on  the  previously  calculated  characteristic,  as  shown 
in  Figure  9. 


!  j 


FIGURE  9.  SHOCK  POINT  CALCULATION 


The  equation  of  the  shock  wave  is  approximated  locally  by: 

-  tan(&!&)  -  ta n(«)  m 

The  properties  at  point  C  are  evaluated  using  the  local  slope  of  the  shock  wave, 
tan  The  point  B  along  the  characteristic  A-D  is  determined  by  the  intersection 
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of  the  6-jj.  characteristic  passing  through  point  A  and  the  6  +  fx  characteristic 
passing  through  point  C.  This  process  is  identical  to  the  field  point  solution  dis¬ 
cussed  previously  except  that  the  properties  at  B  are  determined  by  interpolation 
using  points  E  and  P.  Thus  the  process  is  iterative  in  nature.  Once  the  properties 
at  B  are  determined,  the  compatibility  equation,  Equation  (35),  is  evaluated  along 
B-C  and  the  residual  determined.  Another  shock  angle,  6c,  is  chosen  and  the  pro¬ 
cess  repeated.  Given  the  two  residuals  and  the  two  shock  wave  angles,  an  inter¬ 
polation  can  easily  be  set  up  to  determine  the  shock  angle  for  a  zero  residual. 

Using  this  shock  angle,  the  corresponding  residual  is  calculated  and  these  shock 

angles  and  residuals  are  used  to  generate  yet  another  shock  angle.  This  process 

—6 

continues  until  the  residual  is  less  than  1x10  . 

Since  such  iterative  calculations  are  usually  oscillatory  and  the  first  few  oscil¬ 
lations  may  be  quite  wild,  the  first  shock  angle  used  is  the  shock  angle  that  gener¬ 
ates  sonic  flow  behind  the  shock.  The  second  approximation  is  taken  to  be  the  pre¬ 
vious  shock  angle  &A-  These  two  values  suffice  to  start  the  iteration. 

c.  Body  Point 

The  body  point  solution  musi  satisfy  both  a  geometrical  equation  and  the  con 
patibility  equation.  As  shown  in  Figure  10,  let  point  A  be  a  field  or  shock  point 
and  point  C  be  the  associated  body  point. 


The  line  connecting  A  and  C  is  described  by 


Y-Ya  -  m(X-XA) 


(39) 


For  a  blunt  general  conic  nose,  located  as  shown  in  Figure  11,  the  equation  is: 

Y*  ■  *R(X+a)  -B(X+a)*  (40) 

where  A  is  the  displacement  of  the  nose,  positive  when  the  displacement  is  to  the 
left  of  the  origin. 


FIGURE  11.  BLUNT  BODY  GEOMETRY 


Combining  Equations  (39)  and  (40)  shows  that: 


(41) 


where 

AA  m  B  +  m 

BB  -  2[m(YA-mXA)-7?  +  BA] 

CC  -  (YA-mXA)a  +  a(Ba-2R) 

For  the  case  of  a  body  segment  described  by  a  general  cubic  equation 


Y  -  C,(X-XP)S  4  C*(X-X0)1  +  Cs(X-Xe)  +  c4  (42) 


an  iterative  solution  is  obtained  by  using  the  standard  Newtonian  method: 


xf 


,a-n  .  Y/-"-  (Kf-  Mm  -  Y, 
'c  (dV/tDOc  -  m 


(43) 
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Starting  with  an  approximate  value  of  X^,  values  of  Y^  and  (dY/dX)^  are  deter¬ 
mined,  and  used  in  Equation  (43)  to  obtain  a  better  vaiue  of  X^.  This  process 
continues  until  no  further  significant  change  occurs  in  the  value  of  X^. 

The  compatibility  equation  applied  between  A  and  C  is 

-~2^  (6c  "  *a“  smZsmJx  -  In  ~  «  0  (44) 

For  the  first  approximation,  the  residual  of  Equation  (44)  is  evaluated  as  written 
and  then  Equation  (44)  is  solved  for  P^.  This  new  value  for  now  determines 
a  new  residual.  By  interpolating,  always  using  the  two  latest  new  values  for  Pc 
and  the  associated  residuals,  a  satisfactorily  small  residual  can  be  easily  deter¬ 
mined  in  an  iteration  loop. 

d.  Mesh  Control 

The  method  of  characteristics  program  calculates  along  right  running  charac¬ 
teristics,  i.e. ,  those  that  originate  at  the  shock  wave  and  propagate  downstream 
toward  the  body.  The  X  distance  between  successive  intersections  of  these  char¬ 
acteristics  and  the  body  surface  is  defined  as  the  step  size  of  the  solution.  The 
step  size  is  then  a  function  of  the  length  of  shock  increment  used  to  generate  new 
characteristics.  Subroutine  MESH  both  tests  the  step  size  to  ascertain  its  correct¬ 
ness  and  adjusts  the  shock  wave  step  size  so  that  the  desired  body  increment  is  ob¬ 
tained.  Upon  completion  of  a  characteristic,  MESH  is  used  to  accept  or  reject 
the  characteristic,  depending  on  whether  the  size  criterion  is  satisfied  or  not. 

By  adjusting  the  step  size,  MESH  also  makes  certain  that  there  is  a  characteristic 
that  terminates  at  the  sphere-cone  juncture. 

This  characteristics  program  has  two  distinct  modes  of  operation  in  determin¬ 
ing  the  step  size:  normal  operation  and  juncture  point  operation.  The  program 
switches  automatically  from  one  mode  to  the  other,  as  required.  The  specified 
body  step  size  is  given  the  name  XSTARB.  For  both  the  normal  and  juncture  point 
operation  modes,  MESH  first  calculates  the  step  size,  using  the  body  point  just 
determined,  and  the  body  point  on  the  previously  accepted  characteristic.  This 
candidate  step  size  is  stored  as  DXSTAR,  as  shown  in  Figure  32. 
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FIGURE  12.  SHOCK  AND  BODY  STEP  SIZES 


For  positive  values  of  DXSTAR,  the  program  checks  to  pee  if  DXSTAR  is  more 
than  half  but  less  than  the  full  value  of  XSTARB.  If  DXSTAR  satisfies  these  cri¬ 
teria,  the  complete  characteristic  is  acceptable.  This  is  the  normal  operation. 

For  juncture  point  operation,  the  program  che~k-,  to  see  whether  the  distance  be¬ 
tween  the  b°t  body  point  and  the  juncture  point  is  greater  than  1. 2(XSTARi>;.  If  the 
distance  is  greater,  the  program  returns  control  to  the  main  program,  after  ad¬ 
justing  the  shock  segment  length  to  provide  a  step  size  of  0. 7(XSTARB).  If  the 
distance  is  not  greater  than  1.2 (XSTARB),  the  juncture  point  operation  will  occur 
for  the  next  characteristic.  Initially,  the  shock  wave  segment  length  is  adjusted 
by  multiplying  present  shock  segment  length  by  the  ratio  of  the  desired  body  step 
size  to  that  achieved  using  present  shock  segment  length: 


A®shock, 


ASshocklast 


^  body  juncture  point 
AXl*>dyiast 


On  the  successive  iterative  passes  through  MESH,  which  occur  upon  the  completion 
of  a  new  complete  characteristic,  a  linear  interpolation  is  performed  to  determine 
the  shock  wave  step  size  required  to  place  the  body  point  on  the  juncture  point; 
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shock i  -  shock q]  /  v 

ASshocki+1  »  shock 0  +  AXbodyi  -  AXbodyo  J  V  Abodyjnct pt  "  Xbody0 

(46) 

When  the  program  starts  calculations  from  an  initial  value  line  (IVL),  the  IVL 
itself  sets  its  own  mesh  size.  There  is  no  connection  between  mesh  size  in  the 
IVL  region  and  that  requested  by  the  user  for  the  segment  of  the  body  downstream 
of  the  IVL.  After  the  IVL  has  been  used  to  generate  a  complete  characteristic,  the 
final  step  size  in  the  IVL  region  is  compared  with  that  requested  for  the  appropriate 
body  segment.  If  the  requested  step  size  is  larger  than  that  generated  in  the  IVL 
region,  the  rate  of  growth  of  the  body  step  size  is  limited  to  a  factor  of  2.0. 

e.  Sample  Results 

Figure  13  shows  a  rotationally  symmetric  method  of  characteristics  calculation 
carried  downstream  from  the  IVL  given  in  Figure  7.  Each  point  on  the  IVL  gene¬ 
rates  a  characteristic,  which  moves  inwards  and  downstream  until  it  terminates  on 
the  body.  After  completing  the  IVL-spawned  characteristics,  the  mesh  control 
option  is  brought  into  play.  Under  this  option  the  bow  shock  wave  is  extended  and 
a  characteristic  generated  to  provide  the  desired  step  size  along  the  body. 

These  calculations  are  carried  downstream  until  two  points  on  a  characteristic 
lie  beyond  the  limit  line,  which  is  the  sphere-body  juncture  rotated  clockwise  about 
the  center  of  the  sphere  by  an  angle  equal  to  the  angle  of  attack. 

After  passing  this  limit  line,  another  mode  of  operation  comes  into  play.  This 
mode  takes  the  final  shock  wave  increment  used  to  calculate  a  characteristic  ter¬ 
minating  on  the  body,  and  increases  it  by  10%  for  each  subsequent  characteristic. 
Each  of  these  characteristics  now  terminates  with  a  field  point  when  two  field  points 
lie  downstream  of  the  limit  line.  This  mode  continues  until  a  shock  point  is  calcu¬ 
lated  which  lies  downstream  of  the  limit  line.  This  process  can  be  traced  in 
Figure  13. 


29 


FIGURE  13.  INITIAL  VALUE  SURFACE  FLOW  FIELD 


4.  INITIAL  VALUE  SURFACE  INTERPOLATION  PROGRAM 

Having  calculated  the  complete  flow  field  about  the  spherically  capped,  cone  in  a 
wind-axis  system  an  interpolation  is  carried  out  along  specific  lines  to  determine 
the  properties  in  each  meridian  plane  in  a  body-axis  system.  Each  line  represents 
the  sphere-body  juncture  rotated  clockwise  about  the  center  of  the  sphere  by  an  ( 
angle  /S  which  is  a  function  of  the  desired  meridian  angle  and  the  angle  of  attack. 
When  these  interpolation  lines  are  rotated  about  the  wind  axis,  they  will  form  a 
plane  normal  to  the  body  axis.  ;  ,  ! 

Figure  14  shows  the  constant  body  station  plane  oriented  by  a  wind-axis  system 
at  an  angle  of  attack.  Of .  A  meridian  plane  is  located  by  the  angle  W  in  the  wind- 

i 

axis  system.  From  the  geometry  it  can  be  shown  that: 


COS  &  m  . .  8-.B-.T 

yt+$in*Ytan*  of 

(47) 

sin/3  -  -  sin Y tana, con j8 

(48) 

where  will  vary  between  the,  limits  of  ±01 ;  Note  that  the  angle  £  is  measured 
away  from  the  vertical  axis  in  the  wind -axis  system.  >  1  ! 

»  } 

In  Figure  13,  the  interpolation  lines  are  shown  for  seven  meridian  planes.  ,The 
flow  properties  are  determined  by  linear  interpolation  at  the  intersection  of  each 

t  i 

line  and  each  characteristic.  From  Figure  13  it  is  seen  that  different  lines  end  up 
with  different  numbe.'j  of  interpolated  points.  Using  these  interpolated  data,  a 
second  interpolation  is  carried  but  to  provide  the  same  number  of  points  on  each 
line,  the  number  being  the  number  of  rings  specified'  by  the  user.  The  geometry 
of  the  resulting  surface  for  the  case  shown  in  Figure  13  is  shown  in  Figure  15, 
where  the  number  of  rings  specified  was' five.  This  initial  value  surface  can  be 
used  for  a  blunted  cone  of  25°  half-angle  or  a  slab  delta  wing  with  65°  sweep,  at  a 
Mach  number  of  8. 8  and  an  angle  of  attack  of  15°.  Figure  16  shows ■  anothfer  solu¬ 
tion,  this  time  for  a  blunted,  cone  at  Mach  14.  9  in  helium  at  20°  angle  bf  attack. 

The  interpolation  lines  for  11  meridian  planes  are  shown  in  this  figure,  and  the  re¬ 
sulting  geometry  is  shown  in  Figure  17. 
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SECTION  IV 


THE  THREE-DIMENSIONAL  METHOD  OF  CHARACTERISTICS 

The  original  three-dimensional  method  of  characteristics  has  been  revised  ex¬ 
tensively.  Many  new  approaches  and  procedures  were  used  and  included  (Refs  5,6). 
An  up-to-date  description  of  the  method  and  computer  program  is  presented  here. 

1.  EQUATIONS  AND  ANALYSIS 

The  general  compatibility  relation  for  three-dimensional  steady  flow  of  a  real 
gas  in  equilibrium  is  derived  from  the  fundamental  flow  equations.  In  terms  of 
pressure  and  flow  direction  angles,  this  relation  applies  to  both  isoenergetic  and 
nonisoenergetic  flow. 

a.  Fundamental  Equations 

The  equations  for  steady  inviscid  flow  are 

pV-q  +  q-Vp  =  0  (49) 

pq-Vq  +  VP  -  0  (50) 

q*VS  -  0  (5i) 

where  q  is  the  velocity,  p  the  density,  P  the  pressure  and  s  the  specific  entropy. 
The  gas -dynamic  equation, 

q-7(q2/2)  -  alV  *q  -  0  (52) 

where  a  is  the  speed  of  sound,  is  not  restricted  to  isoenergetic  flow;  however, 
since  it  is  sometimes  derived  for  a  uniform  freestream  (isoenergetic  flow)  only, 
a  general  derivation  is  given  below. 


36 


Consider  P  =  P(p,  8  ) ;  since  VP  =  (0P/0 p)$  Vp  +  (DP/&s)p  Vs  we  have 

q*VP  -  a2qvp  *  0  (53) 

where 

a*  a  CdP/op)s 

by  virtue  of  Equation  (51).  Substitution  of  Equations  (49)  and  (50),  and  the 
identity 

q  vq  -  (vxq)xq  +7Cqa/2) 

into  Equation  (53)  leads  to  the  gas-dynamic  equation  (52).  In  rectangular  coordi¬ 
nates  Equations  (52)  and  (50)  are 

(u^-a^u*  +  (v1-a2)vy  +  (w2-a2)w*  + 

UVCVx+Uy)  +  VWCWy  +  V*)  +  wufu^  +  Wx)  -  o 
UU)!  +  VUy  +  WUj  +  Px/p  -  0 

UVX  +  Wy  +  WVf  +  Py/p  »  0 

UWX  +  VWy  +  WWf  +  Pf/p  -  0 

where  subscripts  x,  y,  z  denote  partial  differentiation. 
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b.  Derivation  of  Compatibility  Relations 

A  characteristic  surface  allows  possible  discontinuities  of  the  first  derivatives 
of  flow  variables  in  the  direction  of  its  normals.  Given  initial  data  on  such  a  sur¬ 
face,  the  normal  derivatives  of  flow  variables  are  not  uniquely  determined  by 
Equation  (54) .  Applying  the  methods  of  Reference  7  by  setting  the  x  and  y  axes 

tangent  to  a  characteristic  surface  at  the  origin  and  solving  for  w  in  terms  of  u, 

z 

v,  w,  P,  and  their  derivatives  with  respect  to  x  and  y,  we  obtain  w  =det[N]/det[D) 

z 

where 
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det  [hi]  m 


U 

1 

0 


v  Gs4-  u4>  ux  +  (a1- vl)  vy  -  uv(v>+uy)  -  w(uw*+vwy) 

0  -UUX  -  VUy  -  ^(/p 

1  “UVX  -  VVy  -  Py/p 


an 


det  [D]  = 


u  v  w2-  a2 

l  o  o 

o  i  o 


Since  the  x  and  y  axes  are  set  tangent  to  a  characteristic  surface,  the  derivative  w 

z 

is  indeterminate.  A  necessary  and  sufficient  condition  for  the  indeterminacy  of  w 

z 

is  det[D]=0  and  det(N]=0.  The  first  leads  to  w  =±a  (the  usual  Mach  cone),  where- 

z 

as  the  second  yields  the  compatibility  relation,  where  w  is  replaced  by  ±a 

a4ux  +  a2vy  *  aCvwy  +  uwx)  +  (vPy  +  uPx)/p  «  0  (55) 

When  the  y-axis  is  chosen  to  coincide  with  a  generatrix  of  the  Mach  cone  (a  bichar¬ 
acteristic)  ,  u  vanishes  and  Equation  (55)  reduces  to 

a4ux  +  a4iy  ¥  avwy  +  vfy/p  -  o 

or 

Ux  *  Vy  ¥  CotMWy  +  3*5^  ?y  -  0  (56) 

pa 


Equation  (56)  is  written  in  a  local,  characteristic-oriented  coordinate  system. 
For  computational  purposes,  it  is  practical  to  relate  it  to  a  fixed  spherical  coordi¬ 
nate  system,  as  shown  in  Figure  18.  The  y-axis  is  replaced  by  the  L-axis  in  a 
bicharacteristic  direction;  the  x  and  z  axes,  by  two  orthogonal  axes  N  and  M,  re¬ 
spectively.  For  a  small,  but  arbitrary,  change  of  velocity  (AC(,  aQ,  aV)  the  cor¬ 
responding  change  of  velocity  components  Au,  Av,  and  Aw  along  the  N,  L,  and 
M  axes  is,  to  the  first  order  of  AC],  and  A"^ : 
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A u  -  qsin6A$  +  qswOcosiAV 

av  s  cos/x a q  -qcosismju.^  +  qsiri&«in66mju 

AW  utsinM&q  *  qco$bcc6ju&&’  T  qsm&sm&casM  AV 

The  partial  derivatives  of  u,  v,  w  with  respect  to  N  or  L  have  the  same  forms  as 
given  above;  hence 

uN  -  qsinfcOM  +  qsin&cos6 I 
VL  -  cos/xqL  -  qeos66in;uOL  +  q6in0sin6sin>u.VL  >  (57) 

WL  »  ±6inju.qu*  qcos6co5>jL0LT  qs»n0sm6cos>jiVL  I 


Substituting  Equation  (57)  into  Equation  (56) ,  we  finally  obtain 


where  subscripts  L,  N  denote  partial  differentiation;  a  is  the  speed  of  sound,  p  the 
density,  jx  the  Mach  angle,  and  6  is  defined  in  Figure  18.  Equation  (49)  is  valid 
for  either  isoenergetic  or  nonisoenergetic  flow;  the  difference  is  the  constants  of 
the  Bernoulli  equation.  For  the  present  application,  air  is  assumed  to  be  a  perfect 
gas  for  which  pa1  is  equal  to  yP ,  where  V  is  the  ratio  of  the  specific  heats. 

Equation  (58)  and  the  Bernoulli  equation  are  used  to  determine  the  supersonic 
flow  field.  A  number  of  numerical  schemes  have  been  proposed  for  solving  Equa¬ 
tion  (58) .  It  can  be  shown  that  many  of  these  schemes  introduce  difficulties  into  the 
solution  by  their  methods  of  approximating  Equation  (58)  in  finite-difference  form 
(Reference  8).  However,  in  this  work  a  simple  form  of  the  generalized  finite- 
difference  approximation  presented  in  Reference  8  is  used. 
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FIGURE  18.  RELATION  BETWEEN  COORDINATE  SYSTEMS 


FIGURE  19.  FIELD  POINT  COMPUTATION 
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c.  Field  Point  Solution 

This  simple  form  of  the  generalized  finite -difference  approximation  is  illustra¬ 
ted  by  the  field  point  solution  shown  in  Figure  19.  A  field  point  is  an  interior  point 
away  from  the  boundaries.  At  such  a  point,  given  an  initial  value  surface,  Equa¬ 
tion  (58)  in  finite -difference  form  along  three  bicharacteristics  passing  through 
this  point  can  be  solved  for  P,  Q  ,  and  V;  however,  in  the  present  formulation, 
one  more  bicharacteristic  is  added  for  accuracy  and  numerical  stability,  as  shown 
in  Figure  19. 

Thus,  in  finite-difference  form,  Equation  (58)  can  be  written 

A'yP0  +  8i&0  +  Ci%  -  Di  {L » 1,4}  (59) 


where 

A{  m  smJlcosjl/lifP) 

Bi  *  **cos6i 
C(  *  sin$6in6i 

Di  *  4iPi  +  Bi«i  '(sm6i&N  +  s«n0cos6iYN)sinjaALi 

and  the  double  summation  convention  is  not  used  in  Dj.  The  barred  quantities  are 
average  properties  between  points  i  and  0,  and  Lj  is  the  distance  between  points 
i  and  0;  6{  is  calculated  using  the  average  flow  direction  between  points  i  and  0 
(see  discussion  in  subsection  2.  c).  As  in  Reference  8  we  define  a  residual  function 

R*  -  E  Uil£  +  Bi©0  DO* 

i»i 

and  seek  a  solution  for  P0,  $o.  and  %  that  minimizes  the  residual  function  ff*. 
Differentiating  U*  with  respect  to  P0,  Bot  and  "%,  and  setting  the  results  equal 
to  zero  yields 
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(60) 


4  4  4  4 

PoZA&i  +«0EBi*  +  %Z*Ci  -  LbiBi 

i«i  i«l  l*i  i«i 


4  4  4  4 

p.EAiQ  +«0EBiCi  +v0ECil  -  EaDi 

l»l  i«i  i«4  i»t 


which  can  then  be  solved  for  P0,  &0<  and 
d.  Body  Point  Solution 

A  body  point  is  a  point  on  the  surface  of  the  fuselage.  At  such  a  point,  as  shown 
in  Figure  20.  the  flow  must  be  tangent  to  the  surface;  hence  at  the  new  body  point  0 

cos  Vq  +  rn^cos&o  +  nbsinB0s\nW0  =  0  (6i) 

where  (1^,  mb>  n^)  are  the  direction  cosines  of  the  normal  vector  to  the  surface 
at  the  new  point.  Since  the  tangency  condition  (61)  relates  the  flow  angles  and 
%,  two  bicharacteristics  are  sufficient  for  determining  P0,  &0,  and  %.  However, 


FIGURE  20.  BODY  POINT  COMPUTATION 


it  is  apparent  from  Figure  20  that  three  bicharacteristics  are  available  and  the 
choice  of  a  specific  pair  will  have  a  strong  influence  on  the  results;  thus  the  addi¬ 
tion  of  a  third  bicharacteristic  is  both  natural  and  desirable. 

This  then  becomes  an  e  itremal  problem  with  the  tangency  condition  (61)  as  a 
constraint.  The  original  procedure  was  to  differentiate  the  residual  function  5* 
with  respect  to  P0  and  considering  ^  as  an  implicit  function  of  %  through 
Equation  (61) ,  and  then  set  the  results  to  zero.  This  led  to  a  transcendental  equa¬ 
tion  which  had  to  be  solved  by  an  iterative  process.  In  the  early  stage  of  revision 
of  the  3DMoC  program,  this  process  occasionally  led  to  an  unwanted  root.  This 
difficulty  was  circumvented  by  a  new  approach  which,  in  effect,  reduces  the  non¬ 
linear  equations  to  linear  ones.  Indeed,  if  we  consider  the  special  case  m^l,  we 
obtain  the  desired  solution  6  =V/ 2  from  Equation  (61)  immediately.  The  finite- 
difference  forms  of  Equation  (58)  :iong  three  bicharacteristics  are  then 

A\ P0  +  Ci%  -  Di-f  Bi  {i*  1.3} 

and  the  minimization  of  the  residual  function 

3 

E*  *  £  (Ail*  4  d-%  -Di  ♦  \  B\ ) 
i«t 

leads  to  the  following  equations  which  are  linear  in  and  % : 


T>ZAil  4 

Y.EAiCi  - 

X>iDi 

iM 

1-t 

i-t 

2  i-t 

3 

3 

3 

T&EAiCi  4 
i-t 

yjtci  = 

i«l 

ECiDi 

a* 

-fEaci 

These  can  readily  be  solved  for  Pe  and  Y£,  For  the  general  case,  m^  =  1,  the  same 
analysis  applies  after  rotating  the  local  coordinates2  to  make  m^l.  After  the  so¬ 
lution  is  obtained  in  terms  of  the  velocity  components  u,  v,  w,  and  the  pressure  P, 
the  solution  is  rotated  back  into  the  original  coordinate  system. 


2  Since  Equation  (50)  relates  the  pressure  P  to  the  local  flow  angles  £  and  V,  it  can 
be  solved  using  any  suitably  chosen  local  coordinates. 
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e.  Shock  Point  Solution 

While  the  field  or  body  point  is  calculated  along  a  streamline,  the  shock  point 
is  calculated  along  a  two-dimensional-shock  line3,  as  shown  in  Figure  21.  At  a 
point  on  the  shock,  the  flow  must  satisfy  the  oblique  shock  relation 


tan  a. 


/  1.-1  \ 

’ 2tfM£>  - - (*+  fit* 

<*+0fo  +  C*-0 

(62) 


where  A0  is  the  flow  deflection  angle  across  the  shock  wave  and  to=Ifc/I£©  is  re¬ 
lated  to  b0  and  V'  through 

is  sin  %  cos  Y0  +  ms  cos  &0  +  ns  sin  &0  sin  %  =  co sAc  (63) 


where  ( 1  ,  m  ,  n  )  are  the  direction  cosines  of  the  freestream  velocity  vector, 
s  s  s 

This  again  is  an  extremal  problem  with  a  constraint  C  (Pe,  43^,  V0)  =0,  given 
parametrically  by  Equations  (62)  and  (63) .  In  the  original  program  the  usual  pro¬ 
cedure  of  minimizing  the  residual  function  led  to  simultaneous  nonlinear  equations 
which  were  difficult  to  solve,  even  numerically.  However,  by  applying  the  same 
technique  of  choosing  appropriate  coordinates,  the  set  of  simultaneous  nonlinear 
equations  can  be  reduced  to  a  single  transcendental  equation,  as  is  shown  below. 

When  the  angle  of  attack  is  zero,  the  freestream  direction  is  aligned  with  the 

Y-axis.  In  this  special  case,  m  =1,  and  $0=  A0is  the  desired  solution  of 

s 

Equation  (63) ;  the  constraint  reduces  to 


tan®°  -  ( y Mif- 1.*  + 1 ) [ 


(*+Oie+(y-i) 


#  Vfz 

5c 


(64) 


relating  fp=  P0/^o  to  %  only.  Three  bicharacteristics  are  used  for  tkie  shock 
point  solution,  as  shown  in  Figure  21,  and  Equation  (58)  assumes  the  finite- 
difference  form 


3A  plane  that  contains  the  freestream  velocity  vector  and  the  downstream  velocity 
vector  at  a  shock  point  is  called  a  "two-dimensional-shock  plane"  (Figure  21b). 
On  the  shock  surface,  a  line  that  is  everywhere  tangent  to  the  two-dimensional- 
shock  planes  is  called  a  "two-dimensional -shock  line". 


A'.P0  +  Bi«0  +  QV.  »  Di 


{i*U} 


where  A^,  B^,  and  are  given  in  Equation  (50) .  Again  we  minimize  the  resi¬ 
dual  function 


3 

EUiP.  +  Bi&o+c^-Di)1 


by  differentiating  it  with  respect  to  P0  and  considering  ^  as  a  function  of  PB 
through  Equation  (64),  and  setting  the  results  to  zero  to  obtain 


E  c Ai  +  Bi  dS/dP)  (AiT> QV. -Di) 
i«4 

tj  £  Ai  q.  +  %  £  Bid  +  P^ECi1  -  £q  q 

t  a  •  A  t _ 4 


-  0 


i»l 


i«t 


i«t 


i«t 


(65) 

(66) 


Equations  (64),  (65),  (66)  are  solved  for  P0,  and  substitution  of  Equations 

(64)  and  (66)  into  (65)  leads  to  a  transcendental  equation  for  1$,  alone,  which  is  then 

solved  by  the  Newton-Raphson  method.  When  the  angle  of  attack  is  not  zero,  the 

coordinates  are  rotated  until  the  local  Y-axis  is  aligned  with  the  freestream  velocity. 

This  makes  m  =  1  and  the  same  analysis  holds.  After  the  solution  has  been  obtained 
s 

in  terms  of  the  velocity  components  u,  v,  w,  and  the  pressure  P,  the  solution  is 
simply  rotated  back  into  the  original  coordinate  system. 

2.  NUMERICAL  PROCEDURES 
a.  Constant  Body  Station  Data  Surfaces 

In  the  original  program,  computation  proceeded  along  outgoing  bicharacteristics 
and  streamlines.  One  drawback  of  such  a  characteristic  network  is  that  the  data 
surfaces  become  more  and  more  distorted  as  the  computation  proceeds  downstream. 
Eventually  the  distortion  (or  warping)  of  the  data  surfaces  makes  their  use  imprac¬ 
tical.  This  happens  much  sooner  in  highly  three-dimensional  flows.  These  warped 
surfaces  also  make  the  interpolation  of  data  or  rearrangement  of  data  points  fairly 
complicated,  as  many  interpolative  operations  in  three  dimensions  must  be  made. 


In  the  new  program,  constant  body  station  data  surfaces  (Y  =  constant)  are  used. 
Each  data  surface  consists  of  a  number  of  data  rings,  with  the  first  data  ring  on  the 
body  and  the  last  data  ring  on  the  shock.  All  data  rings  have  the  same  number  of 
data  points. 

The  advantages  of  the  new  network  are: 

1.  Data  interpolation  is  simple  and  two-dimensional. 

2.  Rearrangement  of  data  points  (which  is  essential  even  in  moderately  three- 
dimensional  flow)  becomes  practical  because  of  simpler  interpolation. 

3.  The  Courant-Friedrichs-Lewy  (CFL)  stability  condition  can  be  readily 
applied  because  the  network  is  two-dimensional  and  more  regular. 

4.  The  extent  of  flow  field  to  be  computed  can  be  easily  determined;  one  does 
not  have  to  figure  out  where  to  begin  the  right  running  characteristics  or 
how  far  a  particular  characteristic  will  reach. 

5.  The  mesh  size  is  more  uniform,  unlike  in  the  old  program  where  the  dis¬ 
tance  between  two  adjacent  data  points  may  become  unusually  large  because 
of  warping  of  the  characteristic  network. 

6.  The  geometric  normal  to  the  shock  surface  is  simple  to  compute  and  this 
makes  the  important  shock  drift  control  practical. 

b.  Base  Point  Location 

Referring  to  Figure  19,  points  1  to  5  are  the  data  points  to  be  used  for  compu¬ 
ting  point  0  on  the  new  data  surface.  Points  1*  to  4'  are  the  base  points,  which  are 
the  intersections  of  the  backward-facing  Mach  conoid  with  the  data  lines.  To  find 
the  location  of  a  base  point  is  one  of  the  basic  operations  of  the  3DMoC  program. 

The  location  of  a  base  point  was  an  involved  and  sometimes  troublesome  opera¬ 
tion  in  the  original  program.  The  coordinates  were  rotated  to  a  local  system  before 
the  intersections  of  the  data  lines  with  the  Mach  cone  were  determined.  The  cor¬ 
rect  choice  of  one  intersection  was  made  and  the  coordinates  were  rotated  back  to 
the  original  system.  No  general  logic  was  available  to  guide  the  choice  of  the  in¬ 
tersection,  which  sometimes  caused  trouble.  The  new  method  described  below 
has  made  possible  a  simple,  direct,  and  unequivocal  determination  of  the  base 
point  location  in  the  new  program. 

Figure  22  shows  data  points  1,  2,  4,  and  5  in  the  upstream  surface,  and  the 
new  point  being  calculated,  denoted  by  0.  The  vector  from  5  to  2  is  "a  and  the 


vector  from  5  to  the  unknown  base  point  2’  is  x .  Point  5  is  the  origin  of  the  stream¬ 
line  through  point  0,  and  is  called  the  hub  point  in  the  following  discussion. 


n 

FIGURE  22.  GEOMETRY  FOR  BASE  POINT  SOLUTION 

Let  x  =  La  ,  where  L  =  x/a  is  an  unknown  to  be  determined.  In  Figure  22  the 
vector  cT  is  defined  by 

C=H-jT  *  k>  -  La  <67) 
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A  =  a2cos2jI  -  (a  - <i|)2 
B  *  a  bcoazjl  -  (a  -q) Cb-q") 
C  =  bacos2M  -  (b-q)2 

The  roots  are 


The  following  is  a  summary  of  the  logic  for  the  choice  of  the  correct  root. 

I.  If  A>0,  the  vector  defined  by  ^"intersects  both  sides  of  the  forward-facing 
Mach  cone.  Two  subcases  must  be  considered. 

(1)  If  C<0,  the  hub  point  {point  5  in  Figure  22)  lies  inside  the  Mach  cone.  The 
positive  root  is  the  desired  one,  and  corresponds  to  the  use  of  the  nega¬ 
tive  sign  in  the  expression  above. 

(2)  If  C  >0,  the  hub  point  lies  outside  the  Mach  cone,  and  two  subcases  are 
possible. 

(a)  If  B<0,  there  are  two  negative  roots.  The  larger  root  is  the  desired 
one,  and  corresponds  to  the  use  of  the  negative  sign. 

(b)  If  B>0,  there  are  two  positive  roots.  Again  the  larger  root  is  the 
desired  one,  and  therefore  the  negative  sign  is  used. 

II.  If  A<0,  the  vector  defined  by  "a" intersects  one  side  of  the  forward -facing  Mach 
cone  and  other  side  of  the  backward -facing  Mach  cone.  Two  subcases  exist. 

(1)  If  C<0,  the  hub  point  lies  inside  Mach  cone,  and  two  subcases  must  be 
considered. 

(a)  If  B<0,  there  are  two  positive  roots.  The  smaller  root  is  wanted. 
Therefore  the  negative  sign  is  used. 

(b)  If  B>0,  there  are  two  negative  roots.  Neither  root  is  acceptable. 

A  failure  message  should  be  printed  and  the  calculation  terminated. 

(2)  If  C>0,  the  hub  point  lies  outside  the  Mach  cone.  There  are  two  sub¬ 
cases. 

(a)  If  B<0,  the  negative  root  is  wanted.  Therefore  the  negative  sign  is 
used. 

(b)  If  3>0,  neither  root  is  acceptable.  A  failure  message  should  be 
printed  and  the  calculation  terminated. 
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ID.  If  A  =  0,  the  vector  defined  by  If  intersects  only  one  side  of  the  forward  Mach 
cone.  The  same  subcases  exist  as  in  n  (i.  e. ,  take  the  negative  sign)  with 
the  exception  that  there  is  only  one  root. 

Therefore,  except  when  A<0  and  B>0  (which  requires  the  selection  of  another  data 
line),  the  negative  sign  is  always  used.  Extrapolation  occurs  when  L  is  not  in  the 
interval  (0,  1).  With  the  factor  L  available,  interpolation  for  daca  at  the  base 
point  is  simple.  For  instance,  the  pressure  P2,  at  point  2'  is  given  by 

P2,-  P5  +  L(Pa-P5)  (72) 

c.  Basic  Solution  Procedure 

The  numerical  solution  is  carried  out  in  a  series  of  data  surfaces  which  are 
perpendicular  to  the  Y-axis.  Given  an  initial  value  surface,  flow  properties  are 
computed  at  points  on  the  next  surface  through  an  iterative  procedure,  following 
streamlines  for  the  body  and  field  points  or  following  two-dimensional -shock  lines 
for  the  shock  points. 

This  procedure  is  illustrated  by  the  field  point  computation.  Referring  to 
Figure  19,  the  location  of  the  new  field  point,  point  0,  is  determined  by  use  of  the 
average  flow  direction  between  point  0  and  the  field  point  on  the  previous  data  sur¬ 
face,  point  5.  In  the  first  iteration,  the  properties  at  point  0  are  assumed  to  be 
the  same  as  those  at  point  5.  Then  the  location  of  the  base  points  1'  to  4*  and  the 
flow  properties  at  these  points  are  determined  by  the  procedure  described  in  the 
preceding  subsection  2.b. 

After  the  four  base  points  are  located,  corresponding  to  assumed  properties 
at  the  new  point  0,  the  coefficients  A.,  B.,  C.,  and  D.  in  Equation  (59)  are  evalu¬ 
ated.  The  angle  6i  (see  Figure  18)  is  determined  by 

cos6i  -  (qixj*N)(|ciixJllN)‘1  (73) 

—  A  -* 

where  N  =  q.  x  L  ;  sin©,  is  given  by 

sin  6i  =  ±  ( l  -  cos26i)  /a  (74) 

which  takes  the  same  sign  as  that  of  X*  (qj  x  f ).  The  derivatives  andT^  are 
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obtained  from 


0N  »  7$  -ft 
*  7V*fi 


The  computation  of  70  and  7V  can  be  illustrated  by  evaluating  70  for  the  first 
bicharacteristic.  We  note  that 


- 

- «»-  ■ 


(75) 


which  are  solved  for  ^ and  the  three  components  of  Vft. 

Oa  aY  at 

Having  evaluated  the  coefficients  A.,  B.,  C.,  and  D.  we  solve  Equation  (60) 
for  the  properties  P0 ,  00 ,  and  which  are  then  used  to  relocate  the  new  field 
point  and  the  corresponding  base  points  for  restarting  the  iteration.  This  process 
is  repeated  until  the  difference  between  two  successive  iterations  is  less  than  a 
specified  value  for  each  of  the  quantities  P0, 0o,  and  Body  points  and  shock 
points  are  computed  in  an  analogous  way. 


d.  Interpolation  of  Data  at  a  Base  Point 

In  the  original  program,  a  second  order  Lagrange  interpolation  was  used  to 
determine  the  flow  properties  at  a  base  point.  At  an  inflection  point  on  a  plot  of 
a  certain  property  versus  distance,  such  an  interpolation  procedure  leads  to  dif¬ 
ferent  values  of  the  property  at  the  same  abscissa,  depending  upon  the  trio  of 
base  points  used,  as  shown  in  Figure  23. 

In  Figure  23,  it  can  be  easily  seen  that  the  different  values  PP.  and  PP0  are 
obtained  by  interpolating,  using  three  points  centered  on  point  in  one  case,  and 
on  point  02  in  the  other.  This  discrepancy  was  responsible  for  some  of  the  troubles 
encountered  in  the  original  program.  In  the  new  program,  linear  interpolation  is 
used  to  determine  the  properties  at  a  base  point.  Comparison  of  numerical  results 
with  experimental  data  has  shown  that  this  simple  interpolation  is  satisfactory. 
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DISTANCE 

FIGURE  23.  INFLECTION  POINT  IN  THE  DATA 

Accuracy  can  be  improved  by  decreasing  the  mesh  size.  Since  a  great  number 
of  iterations  for  properties  at  the  base  points  are  required  in  the  flow  field  compu¬ 
tation,  the  increase  in  computing  time  due  to  the  increased  number  of  data  points 
is  offset  by  the  decrease  in  computing  time  due  to  the  simplicity  of  interpolation. 

e.  beep  Size  Control 

The  Courant-Friedrichs-Lewy  (CFL)  stability  condition,  which  states  that  the 
domain  of  dependence  of  the  difference  equation  must  contain  the  domain  of  depen¬ 
dence  of  the  differential  equation,  is  generally  used  in  numerical  integration  of 
hyperbolic  equations.  Since  the  data  surface  will  be  flat  and  the  network  fairly 
regular,  the  CFL  condition  can  be  applied  easily.  In  the  base  point  computation 
the  position  of  each  base  point  will  be  checked  for  satisfaction  of  the  CFL  condition. 
This  is  simply  done  by  noting  whether  the  factor  L  (see  base  point  computation)  is 
less  than  l/./2«0. 7.  If  this  condition  is  satisfied  at  every  base  point  on  the  data 
surface,  the  step  size  will  be  increased  by  5%;  otherwise,  the  step  size  will  be 
decreased  by  10%. 

The  new  program  also  allows  changes  in  the  grid  size  by  assigning  different 
numbers  of  data  points  in  different  flow  regions.  The  total  number  of  data  points 
per  data  ring  can  be  assigned  for  data  surfaces  between  two  Y  values.  In  addition, 
the  lateral  grid  size  can  be  halved  between  two  assigned  data  points  along  a  data 
ring.  It  was  found  very  useful  to  halve  the  grid  size  in  regions  of  high  lateral 
gradients,  or  high  rates  of  change  of  gradients,  such  as  the  region  near  the  leading 
edge  of  a  wing. 


52 


Notice  that  the  longitudinal  step  size  is  influenced  by  the  lateral  grid  size 
through  the  CFL  condition.  Therefore,  by  assigning  the  number  of  data  points 
and  the  use  of  the  half-size  grid,  the  downstream  step  size  can  be  controlled  in¬ 
directly. 

f.  Respacing  of  Data  Points 

Since  the  network  of  data  points  is  built  up  following  streamlines,  and  stream¬ 
lines  tend  to  converge  on  the  lee  side  of  bodies,  the  data  points  tend  to  crowd  to¬ 
gether  in  the  same  region  after  a  number  of  steps.  Since  the  CFL  condition  has  to 
be  satisfied  everywhere,  this  region  determines  the  step  size.  Even  in  a  moder¬ 
ately  three-dimensional  flow,  the  step  size  may  be  seriously  reduced  after  a  cer¬ 
tain  number  of  steps,  and  a  rearrangement  of  data  points  becomes  necessary.  In 
the  new  program,  the  flat  data  surface  makes  rearranging  data  points  practical. 
Cubic  spline  interpolation  is  used  to  rearrange  the  data  points  twice:  once  along 
the  data  ring  (i. e. ,  in  the  circumferential  direction),  and  once  across  the  data 
rings  (i.e. ,  in  the  radial  direction).  Rearranging  data  points  helps  in  the  numeri¬ 
cal  stability  of  the  computation,  and  also  results  in  a  more  even  distribution  of 
data  points  and  hence  better  accuracy. 

As  the  distance  between  the  shock  and  the  body  increases,  the  distance  between 
two  data  rings  becomes  much  greater  than  the  distance  between  two  data  points. 
When  this  condition  prevails  at  every  shock  point,  the  program  will  automatically 
introduce  an  additional  data  ring  as  the  data  points  are  being  respaced. 

g.  Shock  Point  Drift  Control 

In  three-dimensional  supersonic  inviscid  flow  compuUtions  the  shock  wave  is  a 
free  boundary  to  be  determined  during  the  computation.  As  such,  it  is  subject  to 
drifting.  Experience  has  shown  that,  as  the  number  of  steps  increases,  the  shock 
normal  determined  from  the  local  flow  properties  gradually  deviates  from  the  geo¬ 
metric  shock  normal  until  a  failure  occurs.  This  shock  point  drift  is  shown 
schematically  in  Figure  24. 

In  the  new  program,  subroutine  HARNES  is  used  to  control  the  shock  point 
drift.  After  a  shock  ring  has  been  calculated,  the  shock  normals  determined  from 
the  flow  properties  are  averaged  with  the  geometric  shock  normals  at  every  shock 
point.  The  average  normals  are  then  used  to  recompute  the  flow  properties  at  the 
shock  points.  Experience  has  shown  that,  although  the  average  normal  differs 
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Normal  to  Shock  Line 


FIGURE  24.  SHOCK  POINT  DRIFT 

only  slightly  from  either  of  the  normals  { which  justifies  the  use  of  the  average 
normal) ,  the  shock  point  drift  is  effectively  checked. 

h.  Calculation  of  0N  and  ^  for  the  Field  Point 

In  the  original  program,  a  least  squares  process  was  used  to  obtain  the  average 
gradients  for  £  and  Y ,  using  the  values  of  &  and  Y  at  all  four  base  points  and 
the  new  point.  The  quantities  $N  and  Y#  at  each  base  point  were  then  evaluated 
using  these  average  gradients.  Since,  ideally,  the  local  gradients  of  $  and  Y,  av¬ 
eraged  along  a  given  bicharacteristic,  should  be  used  to  evaluate  and  ^  ,  the 
process  using  only  one  set  of  average  gradients  for  the  entire  Mach  cone  is  inaccu¬ 
rate.  In  the  new  program,  the  local  gradients  of  &  and  Y  are  obtained  using  only 
three  of  the  base  points  and  the  new  point  for  each  bicharacteristic,  as  illustrated 
by  Equation  (75).  Thus  each  bicharacteristic  has  its  own  V&  and  Vf  values  to  be 
used  in  forming  and  Yh. 
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i.  Treatment  of  Undefined  6  and  0 


The  direction  of  the  local  velocity  is  defined  by  the  angles  0  and  ,  while  the 
direction  of  a  bicharacteristic  is  defined  by  the  angle  6  (see  Figure  18) .  The 
angles  b  and  V  become  undefined  when  the  local  velocity  vector  is  aligned  with 
the  Y-axis  (i.e. ,  $■  =0). 

hi  the  original  program  for  computing  the  flow  field  over  a  slab  delta  wing,  this 
occurred  on  the  flat  surface  at  the  plane  of  symmetry  and  was  "fixed"  by  replacing 
the  flat  center  section  by  a  wedge  slope  of  0.01°.  While  this  fix  was  satisfactory 
for  the  slab  delta  wing,  it  was  not  a  general  solution  to  the  problem.  In  the  new 
program,  a  built-in  logic  package  rotates  the  X  and  Y  axes  counter  clockwise 
around  the  Z-axis  through  0. 1  radian  when  the  program  detects  a  small  value  of  & 
The  rotation  could  be  repeated,  when  called  forth  by  a  small  value  of  0,  as  many 
as  16  times.  This  device  has  worked  very  well. 


SECTION  V 


RESULTS 

Both  the  IVS  program  and  the  3DMoC  program  have  been  tested  on  a  number  of 
flow  field  computations.  The  results  of  two  test  cases  are  presented  here.  The 
first  case  is  the  flow  of  helium  over  a  blunted  cone  of  15°  half-angle  at  Mach  14. 9 
and  20°  angle  of  attack.  The  second  case  is  the  flow  field  about  a  slab  delta  wing 
of  70°  sweep  at  Mach  9. 6  and  15°  angle  of  attack.  In  both  cases,  good  agreement 
between  computed  results  and  available  experimental  data  was  obtained. 

1.  BLUNTED  CONE  AT  LARGE  ANGLE  OF  ATTACK 

To  test  the  program's  ability  to  handle  large  angles  of  attack,  the  first  test 
case  is  the  flow  over  a  blunted  cone  of  15°  half-angle  at  20°  angle  of  attack.  The 
freestream  Mach  number  is  14. 9  and  the  medium  is  helium,  which  has  a  ratio  of 
specific  heats  of  1. 6667.  This  is  a  severe  test  case,  since  the  lee  side  of  the  cone 
lies  completely  in  the  shadow  region  of  the  freestream. 

The  body  description  of  a  cone  is  very  simple.  Since  the  cone  is  a  slab  delta 
wing  without  the  flat  sections,  the  same  body  description  subroutine  was  used. 
About  six  minutes  of  CPU  time  were  needed  on  the  CDC  6600  computer  for  calcu¬ 
lating  the  flow  field  to  75  nose  radii  downstream. 

The  computed  shock  shape  is  shown  in  Figure  25,  and  is  compared  with  the 
experimental  data  of  Reference  9.  The  agreement  is  very  good.  Reference  9  also 
presents  the  flow  field  calculations  of  Reference  10  which  failed  at  about  10  nose 
radii  downstream.  Envelope  shocks  may  occur  on  the  lee  side,  but  in  this  case 
they  are  apparently  very  weak  and  the  present  method  is  capable  of  calculating 
right  through  them  without  special  treatment.  Figure  26  is  a  comparison  of  the 
calculated  pressure  with  both  experimental  data  (Reference  9)  and  other  calcu¬ 
lated  data  (Reference  10).  No  essential  disagreement  exist  between  the  two 
methods  up  to  10  nose  radii  downstream,  as  both  agree  quite  well  with  the  experi¬ 
mental  data.  However,  on  the  lee  side,  the  present  method  is  not  restricted  in  its 
downstream  extent. 


FIGURE  25.  SIDE  VIEW  OF  THE  BOW  c,HOCK  ON  A  BLUNTED 
15°  CONE  AT  20°  .  *'  OLE  OF  ATTACK 


I 


1 

2.  SLAB  DELTA  WING  AT  ANGLE  OF  ATTACK 

The  second  test  case  is  the  flow  field  around  a  spherically-capped  70°  sweep 
slab  delta  wing  with  cylindrical  leading  edge's  at. Mach  9.6  and  15°  angle  of  attack. 
The  computation  was  carried  out  to  12.5  nose  radii. downstream,  with  a  CPU  time 
requirement  of  about  three  minutes  on  the  C DC  6600  computer.  Experimental  flata 
were  obtained  from  Bertram's  work  (Reference  11)  for  M  =  9.6  in  air;  Figure  27. 
shows  the  profile  and  the  front  views  of  the  calculated  bow  shock  supported  by  the 
slab  delta  wing.  The  cross  section  is  smooth  and  cbnvex  everywhere.  Figure  28 
shows  the  plan  view  of  the  shock  surface.  Note  that  the  shock  wave  has  a  point  of 

i  • 

inflection  near  the  last  calculated  station.  This  means  that,  somewhat  further 
downstream,  the  bow  shock  will  bend  away  from  the  leading  edge.  This  phenome¬ 
non  can  be  seen  in  Schlieren  photographs  of  Reference  11. 

I 

The  centerline  pressure  distribution  or.  the  windward  surface  is  compared  with 
experimental  data  in  Figure  29.  The  agreement  is  very  good.  Figure  30. compares 

i 

pressure  distributions  along  planes  normal  to  the  wing's  leading  edge  with  experi¬ 
mental  data.  The  abscissa  is  the  ratio  of  distance  along  the  surface  to  wing  thick¬ 
ness.  The  positive  values  of  S/t  refer  to  the  windward  side;  the  negative  values, 
to  the  lee  side.  Very  good  agreement  is  seen,  although  the  Quality  :begins  to  decay 
on  the  expansion  surface  near  the  leading  edge  toward  tHe  end  of  the  wing.  This 
could  be  due  to  viscous  effects.  Nevertheless,  fairly  good  pressure  data  agreement 
is  obtained  on  the  lee  side  of  the  wing  (negative  S/t),  even  thbugh  the  present  numer¬ 
ical  method  does  not  account  for  vortices  which  exist  in  the  real  flow. 
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FIGURE  27.  BOW  SHOCK  ON  A  SPHERICALLY-NOSED  SLAB  DELTA  WING  WITH  70°  SWEEP 


FIGURE  28.  TOP  VIEW  OF  BOW  SHOCK  ON  A  SPHERICALLY-NOSED 

SLAB  DELTA  WING  WITH  70°  SWEEP 


Distance  along  Surface,  S/R 


FIGURE  30.  (Continued) 


FIGURE  30.  (Concluded) 


SECTION  VI 


CONCLUSION 

A  new  version  of  the  three-dimensional  method  cf  characteristics  has  been  de¬ 
veloped  to  calculate  the  steady  supersonic  inviscid  flow  around  smooth  three-dimen¬ 
sional  bodies  at  general  angles  of  attack.  Many  new  approaches  and  procedures 
have  been  devised  and  incorporated  to  improve  the  numerical  method  and  the  com¬ 
puter  program.  Flow  fields  over  various  body  shapes  under  different  freestream 
conditions  have  been  calculated.  The  computed  results  were  compared  with  avail¬ 
able  experimental  data  with  good  agreement  as  exemplified  by  the  two  test  cases 
presented  in  this  report. 

It  has  been  proven  that  the  new  three-dimensional  method  of  characteristics, 
which  is  efficient,  accurate,  and  versatile,  is  an  important  tool  for  flow  field  anal¬ 
ysis  and  vehicle  design.  Possible  extensions  include  the  treatment  of  bodies  with 
surface  slope  discontinuities,  su^n  as  the  compression  or  expansion  comers,  and 
the  calculation  of  embedded  shocks  resulting  from  fins  and  canopies.  The  basic 
schemes  for  calculating  the  embedded  shocks  have  been  conceived  and  are  ready 
for  development. 
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APPENDIX  I 


EVALUATION  OF  THE  COEFFICIENTS  FOR  A  CONIC  SECTION 
Tbe  general  equation  of  a  conic  is 

aY*+  bXY  +  cXa  slY  +  eX  +  f  *  0  <D 

This  can  be  arranged  into  a  more  suitable  form 

Y  «  PX  +  Q  *  ^PX*+SX  +T  (2) 

By  factoring  Equation  (1),  the  following  combination  of  linear  equations  can  be 
derived: 

K(Y-  -KXY-nnaX-ha)*  <Y-m*X-hiMY-m4X-h4)  -  0 


(3) 

where  m^  m2,  — ,  h^,  and  K  . are  constants  to  be  evaluated.  Note  that  Equation  (2) 
has  five  coefficients.  Thus,  a  general  conic  can  be  passed  through  any  five  points 
not  lying  on  a  straight  line;  in  some  cases  the  results  will  be  two  branches  of  a 
hyperbola.  Therefore,  given  five  points  in  a  plane,  as  shown  in  the  sketch  below, 
the  first  four  can  be  connected  by  straight  lines.  These  lines  are  defined  by  the 
slopes ,  m.,  and  the  constants,  h.,  and  these  are  evaluated  using  the  ordinates  and 
abscissas  of  the  four  points. 


Y 


X 


A  ft 


•<D 


Preceding  page  blank 


If  point  1  is  substituted  into  Equation  (3) ,  the  result  is  identically  zero,  since 
point  1  was  used  to  evaluate  the  coefficients  h^,  m 2>  and  hg.  The  same  argu¬ 
ment  holds  for  point  2,  3,  and  4. 

The  value  of  K  can  be  established  by  evaluating  Equation  (3)  at  point  5  to  obtain 
the  following  result 

(Y 5  ~  mzXg  -  ha)(  ^  -  t*n^Xg  -  h4) 

”(YS  -  -  ht>(  Yg  -  msX5  - h5) 

By  expanding  Equation  (3)  and  collecting  like  powers,  it  can  be  shown  that 

P  «  -B/2A 
Q  -  -D/2A 
B  -  ?*-  C/A 
S  =  -PD/A  -  E/A 
T  -  Q*-F/A 

where 

A  «  K  +  l 

B  *  -  [  K  (m4+  +  m4] 

C  «  Km1m5  m2m4 
D  -  -  [ K ( hi  +  Wj)  +  h2  +  h4] 

E  -  K(mth3  +  m3ht)  +  m2h4  +  m4h2 
F  -  Khihj  +  h2h4 

and  the  sign  of  the  radical  must  be  determined  by  evaluation. 

To  determine  the  coefficients  when  four  points  and  one  slope  are  given,  let  the  lo¬ 
cation  of  points  1  and  2  coincide,  with  m g  being  the  given  slope.  The  remainder  of 
the  analysis  follows  identically. 

For  three  points  and  two  slopes,  let  points  1  and  2  coincide  and  points  3  and  4  co¬ 
incide,  with  the  given  slopes  being  mg  and  m^.  The  rest  of  the  analysis  is  the 
same. 


Two  special  forms  are  worth  noting: 


A  straight  line  has  the  coefficients 

P  -  slope 
Q  -  constant 
R  «  S  -  T  -  0 

A  circle  with  its  center  at  the  point  (XQ,  Yq)>  has  the  coefficients 

P-  0 
Q  *  X0 
R  -  -1 
S  -  2Y0 
T  -  P  2  -  Y02 
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LINE  © 

•  Regions  I  and  II :  0  f  Y  <  1 

2*  »  t|-Y4+  2Y 

»  0 

•  Region  m  :  1  s  Y  S  7 

2t  -  1 

Xt  *  (tan2o*)Y  +  (-  tan  20*) 

LINE  © 

•  Region  I:  0  £  Y<  (1  -s*n20°) 

Hz  "  0 

X2  -  }/  -Y*+  2Y 

•  Regions  II  and  HI  :  ( 1  -  sin  20° )  §  Y  £  7 

Z2-  0 

h*  »  (tan20*)Y  i- (cos 20° -tan 20*  + tan 20a«n 20*) 


LINE  © 

•  Regions  I  and  II :  0  5  Y  <  1 

Z,  -  -^-Ya+2Y 
Xs  -  0 

•  Region  III  :  1IYI’ 


lS  -  -1 

X3  -  (tan 20*) v  +  (-tan 20*) 
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Thus,  the  coefficients  P,  Q,  R,  S,  and  T  are  as  follows  : 

j 

•  RegionI:  0  £  Y<  (1  -ein20°) 


Pi  - 

0 

Qt*  0 

a  "i 

st  - 

2 

1 

o 

as,  - 

Pi* 

0 

0 

R2»  0 

S2  * 

0 

Tz-  1 

s<*2  -  o* 

Pa  * 

0 

Qs=  0 

-i 

s3  * 

2 

T*-  0 

SG*  -  -1 

P4  * 

0 

q4*  0 

R4*  0 

s4- 

0 

T4-  i 

SG4  «  o 

Ps- 

0 

QS"  0 

Rs=  0 

s5- 

0 

t5-  1 

1  SG5  -  0 

P*  - 

0 

o 

r6«  0 

s4- 

0 

sg6  -  0 
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•  Region  II:  (1  -sin 20°)  S  Y<  1 

change  to  the  following 

I 

Ps  »  (Um20#)  Q5  *  (cos20°-tan2Cf ♦  fcan20*sn20*) 
Rs»  0  S5»  0  T5*  1  SQ5-  0 

•  Region  HI :  1  S  Y  £  7 

change  to  the  following 


1 

o 

®1*  1 

R,« 

0 

S,  -  0 

T|  -  i  ■ 

ss. 

p3“  0 

Q3=  1 

Rs- 

0 

s5  «  0 

Tj-  1 

SGa 

P4«  (tan2tf) 

Ct4-  (-tan 20 

r4- 

0 

S4-  0 

T*4-  1 

sg4 

(tan  2£f) 

Q4=  C*'  tan  2GT) 

0 

s4-  o 

Tft-  1 

sg4 

4When  the  numerical  value  of  the  square  root  is  identically  zero,  it  is  best  to 
assign  +1  to  T  and  0  to  SG,  so  that  the  SQRT  function  will  not  generate  an 
error  message. 


2.  BLUNTED  15  CONE 


SIDE  VIEW 

Note  that  for  the  blunted  cone,  lines  ©  and  ®  both  lie  in  the  plane  of  symmetry. 

LINE  ® 

i 

•  Region  I : 

0«Y<  (1  -  sin  15°) 

Ht  -  y-Vl  +  ZY 

Xt  =  0 

•  Region  II 

(1  -sinl5°)  5  ys  5 

Z4  *  (tan  15’) Y  +  (cos  15°  -  tan  15*  +  tani5°sinl5°) 

Xt  »  0 

LINE  (J) 

•  Region  I : 

0  *  Y<  (1 -sin  15°) 

•  Region  n  : 


(1-sin  15°)  g  Ys?  5 


Zz  -  o 

X2  -  (tan  15*) V  +  (cos  15° -  tan  15*  +  tanl5°sinl5*) 


LINE  ® 

•  Region  1 :  0  s  Y  <  ( 1  -  sin  15° ) 

Z3  «-• y~Y*  +2Y 
X3  «  0 

•  Region  H:  (1  -sinl5°)  §  Y  £5 

23  -  (-tan  15°) Y  -  (cosl5°- tanl5°  +  tant5°sinl50) 
X3  -  0 


Thus,  the  coefficients  P,  Q,  R,  S,  and  T  are  as  follows: 
•  Region  I:  0  s  Y<  (1 -sinl5°) 


Pi  = 

0 

Qj  *  0 

Ri-  -t 

St  - 

2 

T|«  0 

SGi  a  +1 

Pa  - 

0 

Q2»  0 

0 

S*« 

0 

Tzm  1 

S(a2  *  0 

P5» 

0 

Qs-  C 

P3  =  -1 

S3  = 

2 

Tj®  0 

Sfig  =  ”1 

P4  - 

Q 

Q4  *  0 

r4«  0 

s4  ■ 

0 

T4-  i 

SG4  53  0 

Ps- 

0 

Qg*  0 

Rg»  -1 

Sg* 

2 

t5-  1 

SGg  -  +1 

P*- 

0 

*  0 

R4*  0 

s4* 

0 

T4-  i 

S3*  -  0 

•  Region  II:  (1  -sin  15°)  £  Y  £  5 

change  to  the  following 
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Pi  *  ( tan  15*)  Qim  (cos  15*  -  tan  15°  +  tan  15°stn  15°) 

*t  *  0  5t  -  0  T4  *  1  sat  *  0 


P3-  (-tan  15*) 
Pj  “  o  Ss  - 

P5  -  (tan  15*) 

0  S3  - 


Qs»  -(cos  15*  -tan  15°  +tanl5°s\nl5°) 
o  Ti*  1  SG5  b  0 

Qs  *  (cos  15*  -  tan  15*  +  tani5°sin  15*) 
0  Tg  ■  i  SGs  =*  0 
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j 
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